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Abstract 

We have performed Monte Carlo studies of the 3D random field XY model on L x L x L simple 
cubic lattices, with random field strengths of h r = 1 and 2. We present results for the angle-averaged 
magnetic structure factor, S(k) at L = 64. Our results appear to indicate a phase transition into 
a ferromagnetic state. This is made possible by the existence of a Griffiths singularity. It appears 
that at the phase transition M 2 jumps to zero discontinuously, with a latent heat which is probably 
subextensive. 

PACS numbers: 75.10.Nr, 05.50.+q, 64.60.Cn, 75.10.Hk 
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I. INTRODUCTION 



Some time ago, Monte Carlo calculations 1,3 indicated that there appeared to be a re- 
gion of quasi-long-range order (QLRO) in the phase diagram of the three-dimensional (3D) 
random- field XY model (RFXYM), at weak random field and low temperature. Some recent 
functional renormalization group calculations^ have questioned this. Therefore, since there 
have been substantial improvements in computing power over the last ten years, the author 
felt it worthwhile to conduct a new Monte Carlo study of this model. 

For fixed-length classical spins the Hamiltonian of the RFXYM is 

H = - J^cosOi - - h r ^2 cos(fa - 9i) . (1) 
(y) * 

Each (pi is a dynamical variable which takes on values between and 2tt. The (ij) indicates 

here a sum over nearest neighbors on a simple cubic lattice of size L x L x L. We choose 

each 9i to be an independent identically distributed quenched random variable, with the 

probability distribution 

P(Pi) = 1/27T (2) 

for $i between and 2ir. This Hamiltonian is closely related to models of vortex lattices and 
charge density waves. 

Larking studied a model for a vortex lattice in a superconductor. His model replaces 
the spin-exchange term of the Hamiltonian with a harmonic potential, so that each fa is no 
longer restricted to lie in a compact interval. He argued that for any non-zero value of h r 
this model has no ferromagnetic phase on a lattice whose dimension d is less than or equal 
to four. A more intuitive derivation of this result was given by Imry and Ma,- who assumed 
that the increase in the energy of an L d lattice when the order parameter is twisted at a 
boundary scales as L d ~ 2 . 

As argued by Imry and Ma,- and later justified more carefully,— ^ within a perturbative e 
expansion one finds the phenomenon of "dimensional reduction" . The critical exponents of 
any <i-dimensional O(N) random-field model (where N is the number of components of the 
spin on each site) appear to be identical to those of an ordinary O(N) model of dimension 
d — 2. For the Ising (N = 1) case, this was shown rigorously to be incorrect.- 1 ^ 

Because translation invariance is broken for any non-zero h r , it seems quite implausible 
that the twist energy for Eqn. (1) scales as L d ~ 2 , even though this is correct to all orders in 
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perturbation theory. An alternative derivation by Aizenman and Wehr,— ^ which claims to 
be mathematically rigorous, also makes an assumption equivalent to translation invariance. 
Although the average over the probability distribution of random fields restores translation 
invariance, one must take the infinite volume limit first. It is not correct to interchange the 
infinite volume limit with the average over random fields. 

The basic point is that the Imry-Ma argument for continuous spins (O(N) with N > 2) 
is not self-consistent. One begins by assuming that the random field is weak so that the 
twist energy scales as L d ~ 2 , as in the absence of the random field. Then one shows that if 
d < 4, this cannot be true for large enough L. The only conclusion which should be drawn 
from this is that a deeper analysis of what is going on is needed in that case. 

II. THE GRIFFITHS PHASE 

A mechanism which causes the breakdown of perturbation theory in the presence of 
quenched randomness at the critical temperature, T c , of the pure ferromagnet was discovered 
by Griffiths 1 - in the special case of a bond-diluted Ising model. It was later shown explicitly^ 
that the magnetization per unit volume, M(H) has an essential singularity at H = (where 
H is a uniform field) at and below T c . The result was then extended to the case of binary 
ferromagnetic alloys. 15 It was later argued that this result could be extended beyond the 
Ising model to O(N) models, 16 and that the Griffiths singularity will exist whenever the 
specific heat is divergent at T c for the pure model, so that the bond randomness is a relevant 
perturbation.— 

This argument appears to imply that it is never correct to use a perturbation theory of 
the standard type (e.g. an e expansion) when the randomness is a relevant perturbation. 
On the other hand, an e expansion is only an asymptotic series even in the nonrandom case. 
Thus one may hope that it might continue to be useful even in the presence of the Griffiths 
singularity. 

The original concept of Griffiths^ was a singularity caused by the contributions of large, 
isolated clusters of spins. Hertz, Fleishman and Anderson 1 - developed a more intuitive 
picture of the Griffiths singularity in random bond models, which was developed further 
by Bray and Moore^^. These authors showed that the same type of singular behavior 
would be a consequence of a Lifshits tail of localized eigenvectors of the inverse magnetic 
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susceptibility matrix, near the band edge at eigenvalue zero. 

If one generalizes the idea of the Griffiths singularity from these random bond models to 
random fields models, one finds similarities, but also important differences. The magnetic 
structure factor, S(k) = (|M(k)| 2 ), for XY spins is 

S(k) = £- 3 £cos(k.i^Kcos(&-^)>, ( 3 ) 

i,j 

where r^- is the vector on the lattice which starts at site i and ends at site j, and here the 
angle brackets denote a thermal average. For a random field model, unlike a random bond 
model, the longitudinal part of the magnetic susceptibility, Xi which is given by 

T X (k) = 1-M 2 + L~ 3 ^2 cos(k ■ ?y)((cos(^ — <pj)) - Q {j ) , (4) 

is not the same as S even above T c . For XY spins, 

= (cos(0 i ))(cos(0 j )) + (sin(0 i ))(sin(0 j )) , (5) 

and 

M 2 = L'^Qu = L- 3 ^(cos(0,)) 2 + (sin(0,)> 2 . (6) 

i i 

When there is a ferromagnetic phase transition, S(k = 0) has a stronger divergence than 
x(k = 0). Thus it is S, and not x, which is the best place to look for an essential singularity 
in the random field model. There may also be an effect^ in Xi but we should look in the 
place where the effect is expected to be largest. 

The author is not aware of any explicit studies on the question of the occurrence of 
Lifshits tails in a matrix which has the form of S(k). The natural generalization of the 
Harris criterion is that the Griffiths singularity should be found in S(k = 0) at the T c of the 
pure O(N) ferro magnet for any finite N, and any d > 2. This is because the response of 
the pure system to a weak uniform H field is always divergent at T c , i.e. the susceptibility 
exponent 7 is always positive.— 

A detailed numerical study of the 3D Ising model in a random field (RFIM) at T = 
has been performed by Middleton and Fisher.— These authors find that the 3D RFIM 
has a ferromagnetic critical point, and that the order parameter exponent f3 has a small, 
but positive value. For weak random fields, the 3D RFIM still has two low-temperature 
(ferromagnetic) Gibbs states, even though the random fields have destroyed the Kramers 
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degeneracy. The current author believes the existence of multiple ferromagnetic Gibbs states 
is required for the existence of a critical point with a positive (3. 

Middleton and Fisher find no explicit evidence for the existence of a Griffiths phase in 
the RFIM, although they agree that such a phase should exist. Middleton 23 has also studied 
the 4D RFIM at T = 0. 

The scalar quantity (M 2 ) is a well-defined function of the lattice size L for finite lattices, 
which, with high probability, approaches its large L limit smoothly as L increases. The 
vector M, on the other hand, is not really a well-behaved function of L for an XY model 
in a random field. Knowing the local direction in which M is pointing, averaged over some 
small part of the lattice, may not give us a strong constraint on what (M) for the entire 
lattice will be. 

III. NUMERICAL RESULTS FOR S(k) 

— * 

In this work, we will present results for the average over angles of S'(k), which we write 
as S(k). The data were obtained from L x L x L simple cubic lattices with L = 64 using 
periodic boundary conditions. Some preliminary studies for smaller values of L were also 
done. The calculations were done using a 12-state clock model, i.e. a Z 12 approximation 2 to 
the XY model of Eqn. (1). The strengths of the random field for which data were obtained 
are h r = 1 and 2. These sets of parameters allow the use of a lookup table for the Boltzmann 
factors, because all the energies in the problem are expressible as sums of integers and integer 
multiples of v^3- 24 

This discretization of the phase space of the model has significant effects at very low T, 
but the effects at the temperatures we study here are expected to be negligible compared to 
our statistical errors. The probability distributions for the local magnetization of equilibrium 
states which are calculated for the Z12 model are found to have very small contributions 
from the third and higher harmonics of cos(0) and sin(0). This is strong evidence that the 
12-state clock model is an accurate approximation to the XY model within our range of 
parameters. 

The program uses two linear congruential pseudorandom number generators, one for 
choosing the values of the 9i, and a different one for the Monte Carlo spin flips, which 
are performed by a single-spin-flip heat-bath algorithm. The code was checked by setting 
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FIG. 1: (color online) Angle-averaged structure factor for 64 x 64 x 64 lattices with h r = 1 at 
various temperatures. The error bars indicate one standard deviation statistical errors, and the 
x-axis is scaled logarithmically. 

h r = 0, and seeing that the expected behavior of the pure ferromagnetic system was produced 
correctly. 

Four different L = 64 realizations of the random fields 9i were used. The same four 
samples of random fields were used for all values of T, and both values of h r . Each lattice 
was started off in a random spin state at T/J = 2.3125, significantly above the T c for the pure 
model, and cooled slowly. Thermal averages for S(k) were obtained at a set of temperatures. 
At each T a sequence of 12 spin states obtained at intervals of 20,480 Monte Carlo steps per 
spin (MCS) was Fourier transformed and averaged. The data were then binned according to 
the value of k 2 , to give the angle-averaged S(k). Finally, a logarithically weighted average 
over the four samples was performed. 

The data for h r = 1 at T/J between 1.3125 and 2.0 are shown in Fig. 1. At T/J = 2.0, 
S is clearly flattening out for k < 0.1. As T is lowered, weight is shifted from large values 
of k toward k = 0. For all the smaller values of T/J shown in the figure, the correlation 
length, £, is clearly larger than the size of the samples, L = 64. All of these large £ samples 
show an Ak~ 3 behavior at small k, and the coefficient A is, within our statistical errors, 
independent of T. This apparently universal behavior is consistent with the prediction of 
Giamarchi and Le Doussal 25 for the 3D elastic glass model. For the elastic glass model, this 
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FIG. 2: (color online) Angle-averaged structure factor for 64 x 64 x 64 lattices with h r = 2 at 
various temperatures. The error bars indicate one standard deviation statistical errors, and the 
x-axis is scaled logarithmically. 

Ak~ 3 behavior has be verified numerically by McNamara, Middleton and Zeng.-- For the 
RFXYM, however, we know that this behavior cannot hold down to k = 0, because the 
fixed length of the XY spins yields a sum rule for S: 

£S(k) = L 3 . (7) 

k 

There is also some interesting information to be learned from the dynamical behavior. 
Below T/J= 1.875, one of the four samples shows strong hysteresis, and a second one shows 
mild hysteresis. This may indicate the existence of a first-order phase transition for h r — 1 
at T/J m 1.80. The hysteresis goes away as we continue to still lower T. The value of £ 
under these conditions is clearly bigger than L = 64, however. Therefore we cannot say 
exactly what is going on here from these data. 

Due to the sum rule, Eqn. (7), we know that the Ak~ 3 behavior which is seen in Fig. 1 
cannot be the correct behavior at small k for very large lattices. Increasing the size of our 
lattice by a significant amount is rather impractical. However, it is impossible to distinguish 
between a ferromagnetic state and the QLRO state claimed in the earlier Monte Carlo 
studies^ without studying the true small k limit. 

Another way to learn more about the behavior is to decrease £ by increasing h r . In 
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FIG. 3: (color online) Angle-averaged structure factor for 64 x 64 x 64 lattices with h r = 2 at 
T/J = 0.875, using both hot start and ordered start initial conditions. The error bars indicate one 
standard deviation statistical errors, and the x-axis is scaled logarithmically. 

Fig. 2, we display S(k) data from four h r = 2 samples, obtained using the same procedures 
as before. In this case the peak at k = continues to grow as we lower T, but there does not 
appear to be any simple scaling behavior. The peaks are very flat near k = 0. An attempt 
to fit the data by a Lorentzian shape, or even a Lorentzian raised to some power, is not 
satisfactory. 

At T/J = 0.875, a system which is started in an initially ordered state shows no sign of 
relaxing into the low M 2 state which is found by slow cooling, even on time scales of 10 6 
MCS per spin. In Fig. 3 we show {ln(S(k))), with h r = 2, at this temperature, averaged 
over the same 4 samples for both types of initial conditions. From this figure we see that, 
although the two initial conditions give significantly different values of (M 2 ) the rest of 
(ln(S(k))) is almost the same. The peak in (ln(S(k))) has continued to sharpen slowly with 
the reduction in T. The values of (l^S^O))), not shown in Fig. 3, are 11.54 ± 0.11 for the 
initially ordered state, and 8.26 ± 0.28 for the slow-cooled state. 

The ordered initial conditions give slightly lower values of the energy per spin, as we will 
discuss in more detail in the next section. A careful study to find the highest value of T for 
which the metastability exists was not performed. For h r = 2, unlike h r — 1, the relaxation 
time gets increasingly longer as we continue to lower T beyond the point where metastability 
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first appears. 



IV. DISCUSSION 

If, as argued by Imry and Ma, the random-field energy dominates the spin-exchange 
energy for the 3D RFXYM, how could it be possible to have a ferromagnetic phase? To 
understand this, we need to think carefully about finite-size scaling. 

Consider a very large system at T = 0, with h r — 1 (for example). Since there is no 
degeneracy due to symmetry, and the energy is not quantized, we expect that the ground 
state for any particular sample of random fields should be unique, with probability one. If 
d > 4, then we expect that S'(k) has a 5-function peak at k = 0, of strength M 2 (T = 0) > 0. 
If we consider some low, but finite T, we expect that, in thermodynamic equilibrium, many 
states will contribute to S'(k). These states should have similar values of M 2 > 0, but the 
direction of M ought to be randomly distributed on the unit circle, so that (M) »s 0. 

In d — 3, the Monte Carlo results presented here do not indicate that there are many 
states with similar values of M 2 > 0, but very different values of the direction of M, 
contributing substantially to the equilibrium Gibbs state. This may happen as a rare event 
in a few samples, but it is unlikely to be observed. There is, however, an alternative road 
to ferromagnetism. 

Let us consider some very large, but finite, 3D sample. As Bray has argued,- x _1 nas 
a Lifshits tail of localized states near eigenvalue 0. These localized states occur in regions 
of the sample in which the average local random field has a large magnitude, but a random 
direction. At any temperature, each such localized state corresponds to a compact cluster 
of spins which is pinned to point in a direction which is approximately parallel to the local 
field in that neighborhood. 

When we lower T, the eigenvalue spectrum of x _1 will be modified, as discussed by 
Hertz, Fleishman and Anderson.— For d > 4 mobility edge will move down to zero, and 
the sample polarizes into the renormalized band-edge state. In 3D this type of continuous 
process cannot occur. However, it may become possible, at a low enough T, for the state 
of minimum eigenvalue to polarize most of the sample, without completely destroying the 
other stable clusters. 

A condition which makes it favorable for this process to occur is that the density of 
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localized states should be low, and that the size of a typical localized state should be small 
compared to the average nearest neighbor distance between the localized states. Whether 
this is possible may thus depend on the details of the distribution of random fields. An 
earlier calculation 2 using a different type of random field distribution provided evidence for 
a continuous phase transition. 

A small difference in energy per spin, 

(AE) = (-1.24 ± 0.25) x 1(T 3 J (8) 

per spin, was found in our L = 64 results for h r = 2 at T/J = 0.875 between the slowly 
cooled states and the states relaxed from ordered initial conditions. In order to reach this 
low energy state, the direction of the initial M must be close to the direction of (M) for the 
slow-cooled state. Using an initial M which is not close to this direction results in a state 
which has a higher energy than the slow-cooled state. 

The author believes that this energy is subextensive. This means that, as L — > oo, 
(AE) — ► 0, but L 3 (AE) — > oo. In these terms, the energy difference quoted above is 
—325 ± 66 J. Thus the Boltzmann factor, exp(— 325/0.875), for the weighting of the two 
separate free energy minima is extremely small. Since, of course, the free energies of the 
two Gibbs states must be equal at the phase transition, this Boltzmann factor must be 
compensated by the entropy difference between them. 

It has been known for a long time that something similar occurs 2 ^ 1 ^ in the ID Ising 
model with inverse-square law long-range interactions. For this model it is known that 
the correlation length diverges at T c ,— so that T c is a critical point with /3 — 0. A phase 
transition of this type is also believed to occur in the fc-core percolation model,— which has 
been suggested as model for the ordinary glass transition. 

A rather similar subextensive singularity was recently seen in a type of 2D Ising spin 
glass. 32 If this effect also occurs in the Ising spin glass in higher dimensions, where T c > 0, 
it provides an explanation for the puzzling behavior of the high temperature series for the 
Ising spin glass,— since we know that for random bond Ising models the Griffiths singularity 
appears in 4D. This is consistent with the ideas of Bray and Moore.— ^ 

A phase transition at about T/J« 1.80 for h r — 1 and T j J ~ 0.90 for h r — 2 is in good 
agreement with the work of Gingras and Huse, 1 who worked at T/J = 1.5, and estimated 
the transition to be at h r ~ 1.35. 
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V. SUMMARY 



In this work we have performed Monte Carlo studies of the 3D RFXYM on L = 64 
simple cubic lattices, with random field strengths of h r — 1 and 2. We present results for 
the structure factor, S(k), at a sequence of temperatures. We argue that our results appear 
to indicate a phase transition into a ferromagnetic state. This is made possible by the 
existence of a Griffiths singularity, which invalidates the Imry-Ma analysis. At the phase 
transition M 2 seems to jump to zero discontinuously, but with a latent heat per spin which 
probably goes to zero as L — > oo. 
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